Introduction

This report outlines the initial QC for demodata using the module TAP_H_WTA_v1.0.

For more information about this workflow, please see the GeoMXWorkflows vignette on Bioconductor.

DCCFiles <- list.files(file.path(datadir , DCCdir), pattern=".dcc$", full.names=TRUE)
PKCFiles <- file.path(datadir, PKCfilename)
SampleAnnotationFile <- file.path(datadir, WorkSheet)

myData<-readNanoStringGeoMxSet(dccFiles = DCCFiles,
                                 pkcFiles = PKCFiles,
                                 phenoDataFile = SampleAnnotationFile,
                                 phenoDataSheet = "Template",
                                 phenoDataDccColName = "Sample_ID",
                                 protocolDataColNames = c("aoi", 
                                                          "roi"),
                                 experimentDataColNames = c("panel")) 
## Warning in readNanoStringGeoMxSet(dccFiles = DCCFiles, pkcFiles = PKCFiles, : The following DCC files had no counts: DSP-1001250007864-D-C10.dcc, These will be
## excluded from the GeoMxSet object.
#Shift counts to one to mimic how DSPDA handles zero counts
myData <- shiftCountsOne(myData, elt="exprs", useDALogic=TRUE) 
pkcs <- annotation(myData)
modules <- gsub(".pkc", "", pkcs)

Segment QC

Before excluding any low-performing ROI/AOI segments, we visualize the distributions of the data for the different QC parameters. The cutoffs used are:

Note: You may need to change these cutoffs depending on your experimental setup.

# Default QC cutoffs are commented in () adjacent to the respective parameters
# study-specific values were selected after visualizing the QC results in more
# detail below
QC_params <-
  list(minSegmentReads = 1000, # Minimum number of reads (1000)
       percentTrimmed = 80,    # Minimum % of reads trimmed (80%)
       percentStitched = 80,   # Minimum % of reads stitched (80%)
       percentAligned = 80,    # Minimum % of reads aligned to known targets (80%)
       percentSaturation = 50, # Minimum sequencing saturation (50%)
       minNegativeCount = 10,   # Minimum negative control counts (10)
       maxNTCCount = 1000,     # Maximum counts observed in NTC well (1000)
       minNuclei = 20,         # Minimum # of cells observed in a segment (100)
       minArea = 1000)         # Minimum segment area (5000)
myData <-
  setSegmentQCFlags(myData, 
                    qcCutoffs = QC_params)        

# Collate QC Results
QCResults <- protocolData(myData)[["QCFlags"]]
flag_columns <- colnames(QCResults)
QC_Summary <- data.frame(Pass = colSums(!QCResults[, flag_columns]),
                         Warning = colSums(QCResults[, flag_columns]))
QCResults$QCStatus <- apply(QCResults, 1L, function(x) {
  ifelse(sum(x) == 0L, "PASS", "WARNING")
})
QC_Summary["TOTAL FLAGS", ] <-
  c(sum(QCResults[, "QCStatus"] == "PASS"),
    sum(QCResults[, "QCStatus"] == "WARNING"))
library(ggplot2)

col_by <- "class"

# Graphical summaries of QC statistics plot function
QC_histogram <- function(assay_data = NULL,
                         annotation = NULL,
                         fill_by = NULL,
                         thr = NULL,
                         xlims = NULL) {
  if(is.null(xlims)) {
    xlims <- range(assay_data[[annotation]])
  }
  plt <- ggplot(assay_data,
                aes_string(x = paste0("unlist(`", annotation, "`)"),
                           fill = fill_by)) +
    geom_histogram(bins = 50) +
    geom_vline(xintercept = thr, lty = "dashed", color = "black") +
    theme_bw() + guides(fill = "none") +
    facet_wrap(as.formula(paste("~", fill_by)), nrow = 4) +
    xlim(xlims) +
    labs(x = annotation, y = "Segments, #", title = annotation)
  plt
}

library(knitr)
QC_histogram(sData(myData), "Trimmed (%)", col_by, 80, c(0,101))

QC_histogram(sData(myData), "Stitched (%)", col_by, 80, c(0,101))

QC_histogram(sData(myData), "Aligned (%)", col_by, 75, c(0,101))

QC_histogram(sData(myData), "Saturated (%)", col_by, 50, c(0,101)) +
  labs(title = "Sequencing Saturation (%)",
       x = "Sequencing Saturation (%)")

QC_histogram(sData(myData), "area", col_by, 1000) +
  scale_x_continuous(trans = "log10")

QC_histogram(sData(myData), "nuclei", col_by, 20)

# calculate the negative geometric means for each module
negativeGeoMeans <- 
  esBy(negativeControlSubset(myData), 
       GROUP = "Module", 
       FUN = function(x) { 
         assayDataApply(x, MARGIN = 2, FUN = ngeoMean, elt = "exprs") 
       }) 
protocolData(myData)[["NegGeoMean"]] <- negativeGeoMeans

# explicitly copy the Negative geoMeans from sData to pData
negCols <- paste0("NegGeoMean_", modules)
pData(myData)[, negCols] <- sData(myData)[["NegGeoMean"]]
for(ann in negCols) {
  plt <- QC_histogram(pData(myData), ann, col_by, 2) +
    scale_x_continuous(trans = "log10") 
  print(plt)
}

# detatch neg_geomean columns ahead of aggregateCounts call
pData(myData) <- pData(myData)[, !colnames(pData(myData)) %in% negCols]

Below is a table summarizing the # of segments with a given NTC count:

# show all NTC values, Freq = # of Segments with a given NTC count (for WTA samples this is the deduplicated read count from the empty well A01):
kable(table(NTC_Count = sData(myData)$NTC),
      col.names = c("NTC Count", "# of Segments"), caption = "NTC count summary")
NTC count summary
NTC Count # of Segments
3 64
113 71
397 47
8704 94 
kable(QC_Summary, caption = "QC Summary Table for each Segment")
QC Summary Table for each Segment
Pass Warning
LowReads 272 4
LowTrimmed 276 0
LowStitched 273 3
LowAligned 266 10
LowSaturation 272 4
LowNegatives 66 210
HighNTC 182 94
LowNuclei 276 0
LowArea 265 11
TOTAL FLAGS 42 234

Probe QC

Before we summarize our data into gene-level count data, we will remove low-performing probes. In short, this QC is an outlier removal process, whereby probes are either removed entirely from the study (global) or from specific segments (local). The QC applies to gene targets for which there are multiple distinct probes representing the count for a gene per segment. In WTA data, only one probe exists per target gene; thus, Probe QC does not apply to the endogenous genes in the panel. Rather, it is performed on the negative control probes; there are multiple probes representing our negative controls, which do not target any sequence in the genome. These probes enable calculation of the background per segment and will be important for determining gene detection downstream.

After Probe QC, there will always remain at least one probe representing every gene target. In other words, Probe QC never removes genes from your data.

Set probe QC flags

A probe is removed globally from the dataset if either of the following is true:

  • the geometric mean of that probe’s counts from all segments divided by the geometric mean of all probe counts representing the target from all segments is less than 0.1
  • the probe is an outlier according to the Grubb’s test in at least 20% of the segments

A probe is removed locally (from a given segment) if the probe is an outlier according to the Grubb’s test in that segment.

Note: We do not typically adjust these QC parameters.

# Generally keep the qcCutoffs parameters unchanged. Set removeLocalOutliers to 
# FALSE if you do not want to remove local outliers
myData <- setBioProbeQCFlags(myData, 
                               qcCutoffs = list(minProbeRatio = 0.1,
                                              percentFailGrubbs = 20), 
                               removeLocalOutliers = TRUE)

ProbeQCResults <- fData(myData)[["QCFlags"]]

# Define QC table for Probe QC
qc_df <- data.frame(Passed = sum(rowSums(ProbeQCResults[, -1]) == 0),
                    Global = sum(ProbeQCResults$GlobalGrubbsOutlier),
                    Local = sum(rowSums(ProbeQCResults[, -2:-1]) > 0
                                & !ProbeQCResults$GlobalGrubbsOutlier))

We report the number of global and local outlier probes and exclude those that do not pass Ratio and Global testing.

Probes flagged or passed as outliers
Passed Global Local
18617 1 24 
#Subset object to exclude all that did not pass Ratio & Global testing
cat("Before excluding probes")
## Before excluding probes
dim(myData)
## Features  Samples 
##    18642      276
ProbeQCPassed <- 
  subset(myData, 
         fData(myData)[["QCFlags"]][,c("LowProbeRatio")] == FALSE &
           fData(myData)[["QCFlags"]][,c("GlobalGrubbsOutlier")] == FALSE)

cat("After excluding probes")
## After excluding probes
myData <- ProbeQCPassed 
dim(myData)
## Features  Samples 
##    18641      276

Create gene-level count data

With our Probe QC steps complete, we will generate a gene-level count matrix. The count for any gene with multiple probes per segment is calculated as the geometric mean of those probes. The number below is the number of genes in the raw counts table.

# Check how many unique targets the object has
length(unique(featureData(myData)[["TargetName"]]))
## [1] 18504
# collapse to targets
target_myData <- aggregateCounts(myData)
raw_counts<-data.frame(exprs(target_myData))
raw_counts <- tibble::rownames_to_column(raw_counts, "Gene")
out_raw_counts<-paste0(output_prefix,'_raw_counts.csv')
write.table(raw_counts,out_raw_counts,sep=",",row.names = F,col.names=T,quote=F)

Segment and gene filtering

In addition to Segment and Probe QC, we also determine the limit of quantification (LOQ) per segment. The LOQ is calculated based on the distribution of negative control probes and is intended to approximate the quantifiable limit of gene expression per segment. Please note that this process is more stable in larger segments. Likewise, the LOQ may not be as accurately reflective of true signal detection rates in segments with low negative probes counts (ex: <2). The formula for calculating the LOQ in a \(i^{th}\) segment is:

\[LOQ_{i} = geomean(NegProbe_{i}) * geoSD(NegProbe_{i})^{n}\]

Nanostring typically uses 2 geometric standard deviations (\(n = 2\)) above the geometric mean as the LOQ, which is reasonable for most studies, and recommend that a minimum LOQ of 2 be used if the LOQ calculated in a segment is below this threshold.

After determining the limit of quantification (LOQ) per segment, we recommend filtering out either segments and/or genes with abnormally low signal (below the LOQ). Filtering is an important step to focus on the true biological data of interest.

# Define LOQ SD threshold and minimum value
cutoff <- 2
minLOQ <- 2

# Calculate LOQ per module tested
LOQ <- data.frame(row.names = colnames(target_myData))
for(module in modules) {
  vars <- paste0(c("NegGeoMean_", "NegGeoSD_"),
                 module)
  if(all(vars[1:2] %in% colnames(pData(target_myData)))) {
    LOQ[, module] <-
      pmax(minLOQ,
           pData(target_myData)[, vars[1]] * 
             pData(target_myData)[, vars[2]] ^ cutoff)
  }
}
pData(target_myData)$LOQ <- LOQ

LOQ_Mat <- c()
for(module in modules) {
  ind <- fData(target_myData)$Module == module
  Mat_i <- t(esApply(target_myData[ind, ], MARGIN = 1,
                   FUN = function(x) {
                     x > LOQ[, module]
                   }))
  LOQ_Mat <- rbind(LOQ_Mat, Mat_i)
}
# ensure ordering since this is stored outside of the geomxSet
LOQ_Mat <- LOQ_Mat[fData(target_myData)$TargetName, ]

Segment filtering

We first filter out segments with exceptionally low signal. These segments will have a small fraction of panel genes detected above the LOQ relative to the other segments in the study.

# Save detection rate information to pheno data
pData(target_myData)$GenesDetected <- 
  colSums(LOQ_Mat, na.rm = TRUE)
pData(target_myData)$GeneDetectionRate <-
  pData(target_myData)$GenesDetected / nrow(target_myData)

# Determine detection thresholds: 1%, 5%, 10%, 15%, >15%
pData(target_myData)$DetectionThreshold <- 
  cut(pData(target_myData)$GeneDetectionRate,
      breaks = c(0, 0.01, 0.05, 0.1, 0.15, 1),
      labels = c("<1%", "1-5%", "5-10%", "10-15%", ">15%"))

# stacked barplot of different cutpoints (1%, 5%, 10%, 15%)
ggplot(pData(target_myData),
       aes(x = DetectionThreshold)) +
  geom_bar(aes(fill = `class`)) +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
  theme_bw() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(x = "Gene Detection Rate",
       y = "Segments, #",
       fill = "Segment class")

Table showing the tissue classes impacted by each threshold:

# cut percent genes detected at 1, 5, 10, 15
kable(table(pData(target_myData)$DetectionThreshold,
            pData(target_myData)$class))
DKD normal
<1% 7 1
1-5% 3 6
5-10% 15 7
10-15% 28 5
>15% 106 94

Generally, 5-10% detection is a reasonable segment filtering threshold. For now, we will select segments that have a 5% detection rate.

Note: These guidelines may require adjustment depending on the experimental design (e.g. segment types, size, nuclei) and tissue characteristics (e.g. type, age).

geneDetectionRateThresh<-0.05

target_myData <-
  target_myData[, pData(target_myData)$GeneDetectionRate >= geneDetectionRateThresh]

dim(target_myData)
## Features  Samples 
##    18504      255
library(scales) # for precent

# Calculate detection rate:
LOQ_Mat <- LOQ_Mat[, colnames(target_myData)]
fData(target_myData)$DetectedSegments <- rowSums(LOQ_Mat, na.rm = TRUE)
fData(target_myData)$DetectionRate <-
  fData(target_myData)$DetectedSegments / nrow(pData(target_myData))

Gene filtering

We will graph the total number of genes detected in different percentages of segments. Based on the visualization below, we can better understand global gene detection in our study and select how many low detected genes to filter out of the dataset. Gene filtering increases performance of downstream statistical tests and improves interpretation of true biological signal.

# Plot detection rate:
plot_detect <- data.frame(Freq = c(1, 5, 10, 20, 30, 50))
plot_detect$Number <-
  unlist(lapply(c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5),
                function(x) {sum(fData(target_myData)$DetectionRate >= x)}))
plot_detect$Rate <- plot_detect$Number / nrow(fData(target_myData))
rownames(plot_detect) <- plot_detect$Freq

ggplot(plot_detect, aes(x = as.factor(Freq), y = Rate, fill = Rate)) +
    geom_bar(stat = "identity") +
    geom_text(aes(label = formatC(Number, format = "d", big.mark = ",")),
              vjust = 1.6, color = "black", size = 4) +
    scale_fill_gradient2(low = "orange2", mid = "lightblue",
                         high = "dodgerblue3", midpoint = 0.65,
                         limits = c(0,1),
                         labels = scales::percent) +
    theme_bw() +
    scale_y_continuous(labels = scales::percent, limits = c(0,1),
                       expand = expansion(mult = c(0, 0))) +
    labs(x = "% of Segments",
         y = "Genes Detected, % of Panel > LOQ")

We will now select genes detected in at least 10% of segments and filter out the remainder of the targets.

Note: if we know that a key gene is represented in only a small number of segments (<10%) due to biological diversity, we may select a different cutoff or keep the target gene by manually selecting them for inclusion in the data object.

# Subset to target genes detected in at least 1 of the samples.
#   Also manually include the negative control probe, for downstream use
targetDetectionRateThresh<-0.10
negativeProbefData <- subset(fData(target_myData), CodeClass == "Negative")
neg_probes <- unique(negativeProbefData$TargetName)
target_myData <- 
   target_myData[fData(target_myData)$DetectionRate > targetDetectionRateThresh |
                     fData(target_myData)$TargetName %in% neg_probes, ]

"Remaining targets and ROIs"
## [1] "Remaining targets and ROIs"
dim(target_myData)
## Features  Samples 
##     9906      255

Normalization

We will now normalize the GeoMx data for downstream visualizations and differential expression. The two common methods for normalization of DSP-NGS RNA data are i) quartile 3 (Q3) or ii) background normalization.

Q3 normalization is typically the preferred normalization strategy for most DSP-NGS RNA studies, so we will perform Q3 normalization here.

Before normalization, we will explore the relationship between the upper quartile (Q3) of the counts in each segment with the geometric mean of the negative control probes in the data. Ideally, there should be a separation between these two values to ensure we have stable measure of Q3 signal. If you do not see sufficient separation between these values, you may consider more aggressive filtering of low signal segments/genes. For additional conceptual guidance, please refer to Nanostring’s Data Analysis White Paper for DSP-NGS Assays.

library(reshape2)  # for melt
library(cowplot)   # for plot_grid

# Graph Q3 value vs negGeoMean of Negatives
ann_of_interest <- "class"
Stat_data <- 
  data.frame(row.names = colnames(exprs(target_myData)),
             Segment = colnames(exprs(target_myData)),
             Annotation = pData(target_myData)[, ann_of_interest],
             Q3 = unlist(apply(exprs(target_myData), 2,
                               quantile, 0.75, na.rm = TRUE)),
             NegProbe = exprs(target_myData)[neg_probes, ])
Stat_data_m <- melt(Stat_data, measure.vars = c("Q3", "NegProbe"),
                    variable.name = "Statistic", value.name = "Value")

plt1 <- ggplot(Stat_data_m,
               aes(x = Value, fill = Statistic)) +
  geom_histogram(bins = 40) + theme_bw() +
  scale_x_continuous(trans = "log2") +
  facet_wrap(~Annotation, ncol = 3) + 
  scale_fill_brewer(palette = 3, type = "qual") +
  labs(x = "Counts", y = "Segments, #")

plt2 <- ggplot(Stat_data,
               aes(x = NegProbe, y = Q3, color = Annotation)) +
  geom_abline(intercept = 0, slope = 1, lty = "dashed", color = "darkgray") +
  geom_point() + guides(color = "none") + theme_bw() +
  scale_x_continuous(trans = "log2") + 
  scale_y_continuous(trans = "log2") +
  theme(aspect.ratio = 1) +
  labs(x = "Negative Probe GeoMean, Counts", y = "Q3 Value, Counts")

plt3 <- ggplot(Stat_data,
               aes(x = NegProbe, y = Q3 / NegProbe, color = Annotation)) +
  geom_hline(yintercept = 1, lty = "dashed", color = "darkgray") +
  geom_point() + theme_bw() +
  scale_x_continuous(trans = "log2") + 
  scale_y_continuous(trans = "log2") +
  theme(aspect.ratio = 1) +
  labs(x = "Negative Probe GeoMean, Counts", y = "Q3/NegProbe Value, Counts")

btm_row <- plot_grid(plt2, plt3, nrow = 1, labels = c("B", ""),
                     rel_widths = c(0.43,0.57))
plot_grid(plt1, btm_row, ncol = 1, labels = c("A", ""))

# Q3 norm (75th percentile) for WTA/CTA  with or without custom spike-ins
target_myData <- normalize(target_myData , data_type = "RNA",
                             norm_method = "quant", 
                             desiredQuantile = .75,
                             toElt = "q_norm")

## save q3 data as csv file

norm_data<-data.frame(assayDataElement(target_myData , elt = "q_norm"))
norm_data <- tibble::rownames_to_column(norm_data, "Gene")
out_norm_data<-paste0(output_prefix,'_Q3_normdata.',
                      'geneDetect',
                      geneDetectionRateThresh,
                      '_',
                      'targetDetect',
                      targetDetectionRateThresh,
                      '.csv')

write.table(norm_data,out_norm_data,sep=",",row.names = F,col.names=T,quote=F)
# lets also save the entire targetmyData as an rstructure
save(target_myData, file = paste0(output_prefix,".RData"))

To demonstrate the effects of normalization, we graph representative boxplots of the data for individual segments before and after normalization.

# visualize the first 10 segments with each normalization method
boxplot(exprs(target_myData)[,1:20],
        col = "#9EDAE5", main = "Raw Counts",
        log = "y", names = 1:20, xlab = "Segment",
        ylab = "Counts, Raw",
        ylim=c(1,40000))

boxplot(assayDataElement(target_myData[,1:20], elt = "q_norm"),
        col = "#2CA02C", main = "Q3 Norm Counts",
        log = "y", names = 1:20, xlab = "Segment",
        ylab = "Counts, Q3 Normalized",
         ylim=c(1,15000))

Unsupervised analysis

Now that we have filtered and normalized our data set we can explore the data.

UMAP & t-SNE & PCA

library(umap)
library(Rtsne)

myshapes<-c(16,17,18,15,21,22,3,42,4,8)
# update defaults for umap to contain a stable random_state (seed)
custom_umap <- umap::umap.defaults
custom_umap$random_state <- 42
# run UMAP
umap_out <-
  umap(t(log2(assayDataElement(target_myData , elt = "q_norm"))),  
       config = custom_umap)
pData(target_myData)[, c("UMAP1", "UMAP2")] <- umap_out$layout[, c(1,2)]
ggplot(pData(target_myData),
       aes(x = UMAP1, y = UMAP2, color = class, shape=region)) +
  geom_point(size = 3) +
  scale_shape_manual(values=myshapes) +
  theme_bw()

# run tSNE
set.seed(42) # set the seed for tSNE as well
tsne_out <-
  Rtsne(t(log2(assayDataElement(target_myData , elt = "q_norm"))),
        perplexity = ncol(target_myData)*.15)
pData(target_myData)[, c("tSNE1", "tSNE2")] <- tsne_out$Y[, c(1,2)]
ggplot(pData(target_myData),
       aes(x = tSNE1, y = tSNE2, color = class, shape=region)) +
  geom_point(size = 3) +
  scale_shape_manual(values=myshapes) +
  theme_bw()

## run PCA
PCAx<-1
PCAy<-2
PCAxy <- c(as.integer( PCAx ),as.integer( PCAy) ) # selected principal components


pca.object <- prcomp(t(log2(assayDataElement(target_myData , elt = "q_norm"))))
pcaData = as.data.frame(pca.object$x[, PCAxy]); 
pData(target_myData)[, c("PC1", "PC2")] <- pcaData[,c(1,2)]
percentVar=round(100*summary(pca.object)$importance[2, PCAxy],0)


ggplot(pData(target_myData),
       aes(x = PC1, y = PC2, color = class, shape=region)) +
  geom_point(size = 3) +
  xlab(paste0("PC", PCAx ,": ", percentVar[1], "% variance")) +
  ylab(paste0("PC", PCAy ,": ", percentVar[2], "% variance")) +
  scale_shape_manual(values=myshapes) +
  theme_bw()

Heat map of high CV (>0.8) genes

library(pheatmap)  # for pheatmap
# create a log2 transform of the data for analysis
assayDataElement(object = target_myData, elt = "log_q") <-
  assayDataApply(target_myData, 2, FUN = log, base = 2, elt = "q_norm")

# create CV function
calc_CV <- function(x) {sd(x) / mean(x)}
CV_dat <- assayDataApply(target_myData,
                         elt = "log_q", MARGIN = 1, calc_CV)

# Identify genes in the top 3rd of the CV values
GOI <- names(CV_dat)[CV_dat > quantile(CV_dat, 0.8)]


pheatmap(assayDataElement(target_myData[GOI, ], elt = "log_q"),
         scale = "row", 
         show_rownames = FALSE, show_colnames = FALSE,
         border_color = NA,
         clustering_method = "average",
         clustering_distance_rows = "correlation",
         clustering_distance_cols = "correlation",
         breaks = seq(-3, 3, 0.05),
         color = colorRampPalette(c("purple3", "black", "yellow2"))(120),
         annotation_col = pData(target_myData)[, c("class", "region")])

R tool version information

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] pheatmap_1.0.12         Rtsne_0.16              umap_0.2.8.0            cowplot_1.1.1           reshape2_1.4.4          scales_1.2.0           
##  [7] knitr_1.38              GeomxTools_1.99.4       NanoStringNCTools_1.1.2 S4Vectors_0.30.2        Biobase_2.52.0          BiocGenerics_0.38.0    
## [13] RColorBrewer_1.1-3      ggplot2_3.3.5           colortools_0.1.5       
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.0           uuid_1.1-0             shadowtext_0.1.1       fastmatch_1.1-3        systemfonts_1.0.4      plyr_1.8.7            
##   [7] igraph_1.3.0           lazyeval_0.2.2         splines_4.1.1          BiocParallel_1.26.2    GenomeInfoDb_1.28.4    digest_0.6.29         
##  [13] yulab.utils_0.0.4      htmltools_0.5.2        GOSemSim_2.18.1        viridis_0.6.2          GO.db_3.13.0           lmerTest_3.1-3        
##  [19] fansi_1.0.3            magrittr_2.0.3         memoise_2.0.1          cluster_2.1.3          Biostrings_2.60.2      graphlayouts_0.8.0    
##  [25] askpass_1.1            enrichplot_1.12.3      colorspace_2.0-3       blob_1.2.3             ggrepel_0.9.1          textshaping_0.3.6     
##  [31] xfun_0.30              dplyr_1.0.8            EnvStats_2.7.0         crayon_1.5.1           RCurl_1.98-1.6         jsonlite_1.8.0        
##  [37] scatterpie_0.1.7       lme4_1.1-29            ape_5.6-2              glue_1.6.2             polyclip_1.10-0        gtable_0.3.0          
##  [43] zlibbioc_1.38.0        XVector_0.32.0         V8_4.1.0               DOSE_3.18.3            DBI_1.1.2              randomcoloR_1.1.0.1   
##  [49] ggthemes_4.2.4         Rcpp_1.0.8.3           viridisLite_0.4.0      reticulate_1.24        gridGraphics_0.5-1     tidytree_0.3.9        
##  [55] bit_4.0.4              htmlwidgets_1.5.4      httr_1.4.2             fgsea_1.18.0           ellipsis_0.3.2         pkgconfig_2.0.3       
##  [61] farver_2.1.0           sass_0.4.1             utf8_1.2.2             ggplotify_0.1.0        tidyselect_1.1.2       labeling_0.4.2        
##  [67] rlang_1.0.2            AnnotationDbi_1.54.1   cellranger_1.1.0       munsell_0.5.0          tools_4.1.1            cachem_1.0.6          
##  [73] downloader_0.4         cli_3.2.0              generics_0.1.2         RSQLite_2.2.12         evaluate_0.15          stringr_1.4.0         
##  [79] fastmap_1.1.0          yaml_2.3.5             ragg_1.2.2             ggtree_3.0.4           outliers_0.15          bit64_4.0.5           
##  [85] tidygraph_1.2.1        purrr_0.3.4            KEGGREST_1.32.0        ggraph_2.0.5           nlme_3.1-157           ggiraph_0.8.2         
##  [91] aplot_0.1.3            DO.db_2.9              compiler_4.1.1         rstudioapi_0.13        beeswarm_0.4.0         curl_4.3.2            
##  [97] png_0.1-7              treeio_1.16.2          tibble_3.1.6           tweenr_1.0.2           bslib_0.3.1            stringi_1.7.6         
## [103] highr_0.9              RSpectra_0.16-0        lattice_0.20-45        Matrix_1.4-1           nloptr_2.0.0           vctrs_0.4.1           
## [109] pillar_1.7.0           lifecycle_1.0.1        jquerylib_0.1.4        data.table_1.14.2      bitops_1.0-7           patchwork_1.1.1       
## [115] qvalue_2.24.0          R6_2.5.1               gridExtra_2.3          vipor_0.4.5            IRanges_2.26.0         boot_1.3-28           
## [121] MASS_7.3-56            assertthat_0.2.1       openssl_2.0.0          rjson_0.2.21           withr_2.5.0            GenomeInfoDbData_1.2.6
## [127] clusterProfiler_4.0.5  grid_4.1.1             ggfun_0.0.6            minqa_1.2.4            tidyr_1.2.0            rmarkdown_2.13        
## [133] ggforce_0.3.3          numDeriv_2016.8-1.1    tinytex_0.38           ggbeeswarm_0.6.0